//
//  Jacobian.cpp
//  EMT_C
//
//  Created by Duong Ngo on 8/27/16.
//  Copyright © 2016 Duong Ngo. All rights reserved.
//

#include "Jacobian.hpp"
#include "mapping.hpp"

#include <math.h>

void find_jacobian_structure_for_time_t(const int& t, Index* iRow, Index* jCol, int& dem)
{
    // Mapping f
    std::map<std::string, int> f;
    find_map_at_time_t(t, f);
    
    
    
    //Bankers
    //g0
    int pg= (t-1)*var_foc;
    iRow[dem]=pg; jCol[dem]=f.at("s.gamma"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.cb"); ++dem;
    
    //g1
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.gamma"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.Rf"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.gamma_f"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.ip_f"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.muc"); ++dem;
    
    //g2
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.gamma"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.Rm"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.gamma_f"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.ip_f"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.muc"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.mur"); ++dem;
    
    //g3
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.gamma"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.gamma_f"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.ip_f"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.muc"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.mur"); ++dem;

    //g4
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.gamma"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.ql"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.ql_f"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.gamma_f"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.ip_f"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.muc"); ++dem;

    //g5
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.n_p"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.ip"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.tau"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.n"); ++dem;
    
    

    //g6
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.m"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.Rm_p"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.m_p"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.ip"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.ql"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.s"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.bh"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.bh_p"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.cb"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.tau"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.n_p"); ++dem;

    //g7
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.n"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.m"); ++dem;
    
    //g8
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.n"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.bh"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.m"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.muc"); ++dem;
    
    //g9
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.bh"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.bh_p"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.ip"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.s"); ++dem;

   
    //Households
    
    //g10
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.lama"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.etaz"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.ch"); ++dem;
 
    //g11
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.lamb"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.ch"); ++dem;
    
    //g12
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.lama"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.Rm"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.lamb_f"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.ip_f"); ++dem;

    //g13
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.ql"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.lamb"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.ql_f"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.lamb_f"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.ip_f"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.etab"); ++dem;

    //g14
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.lamb"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.lamb_f"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.pm_f"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.lama_f"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.k"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.y_f"); ++dem;

    //g15
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.Rm_p"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.m_p"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.ip"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.ql"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.s"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.tau"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.ch"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.i"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.bh_p"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.etaz"); ++dem;
    
    
    //g16
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.bh"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.etab"); ++dem;
    
    //g17
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.pie"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.lama_f"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.lama"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.pie_f"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.y_f"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.y"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.pm"); ++dem;

    //g18
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.y"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.k_p"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.l"); ++dem;

    
    //g19
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.u"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.pie_f"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.y"); ++dem;

  
    //g20
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.y"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.cb"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.ch"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.i"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.bh"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.pie"); ++dem;
    
    //g21
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.k"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.k_p"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.i"); ++dem;
    
    
    //g22
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.ip"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.pie"); ++dem;

    //g23
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.Rf"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.u"); ++dem;
    
    //g24
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.l"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.pm"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.y"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.lama"); ++dem;

}

//***********************************************************************************
//***********************************************************************************

void find_jacobian_values_for_time_t (const int& t, const std::vector<double>& ss, const Number* x, Number* values, int& dem)
{
    // Var s
    variables s;
    find_variable_s_at_time_t(t, x, s);
    
    //Derivatives
    //Bankers
    //g0
    int pg= (t-1)*var_foc;
//    iRow[dem]=pg; jCol[dem]=f.at("s.gamma"); ++dem;
    values[dem]= s.cb;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.cb"); ++dem;
    values[dem]= s.gamma;
    ++dem;
    
    //g1
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.gamma"); ++dem;
    values[dem]= 1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.Rf"); ++dem;
    values[dem]= -beta_b*s.gamma_f*s.ip_f;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.gamma_f"); ++dem;
    values[dem]= -beta_b*s.Rf*s.ip_f;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.ip_f"); ++dem;
    values[dem]= -beta_b*s.Rf*s.gamma_f;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.muc"); ++dem;
    values[dem]= -2*pe*pow(fmax(s.muc,0),1);
    ++dem;
    
    //g2
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.gamma"); ++dem;
    values[dem]= 1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.Rm"); ++dem;
    values[dem]= -beta_b*s.gamma_f*s.ip_f;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.gamma_f"); ++dem;
    values[dem]= -beta_b*s.Rm*s.ip_f;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.ip_f"); ++dem;
    values[dem]= -beta_b*s.Rm*s.gamma_f;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.muc"); ++dem;
    values[dem]= -2*pe*pow(fmax(s.muc,0),1);
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.mur"); ++dem;
    values[dem]= -varphi;
    ++dem;
    
    //g3
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.gamma"); ++dem;
    values[dem]= 1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.gamma_f"); ++dem;
    values[dem]= -beta_b*Rn*s.ip_f;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.ip_f"); ++dem;
    values[dem]= -beta_b*Rn*s.gamma_f;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.muc"); ++dem;
    values[dem]= -2*pe*pow(fmax(s.muc,0),1);
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.mur"); ++dem;
    values[dem]= -1;
    ++dem;
    
    //g4
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.gamma"); ++dem;
    values[dem]= s.ql+thetaa;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.ql"); ++dem;
    values[dem]= s.gamma;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.ql_f"); ++dem;
    values[dem]= -beta_b*deltab*s.gamma_f*s.ip_f;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.gamma_f"); ++dem;
    values[dem]= -beta_b*(deltab+ deltab*s.ql_f)*s.ip_f;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.ip_f"); ++dem;
    values[dem]= -beta_b*(deltab+ deltab*s.ql_f)*s.gamma_f;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.muc"); ++dem;
    values[dem]= -2*(1-vec_kappa[t])*pe*pow(fmax(s.muc,0),1);
    ++dem;
    
    //g5
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.n_p"); ++dem;
    values[dem]= s.ip;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.ip"); ++dem;
    values[dem]= s.n_p;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.tau"); ++dem;
    values[dem]= 1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.n"); ++dem;
    values[dem]= -1;
    ++dem;
    
    //g6
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.m"); ++dem;
    values[dem]= 1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.Rm_p"); ++dem;
    values[dem]= -s.m_p*s.ip;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.m_p"); ++dem;
    values[dem]= -s.Rm_p*s.ip;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.ip"); ++dem;
    values[dem]= -s.Rm_p*s.m_p+ deltab*s.bh_p+ (Rn-1)*s.n_p;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.ql"); ++dem;
    values[dem]= -s.s;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.s"); ++dem;
    values[dem]= -s.ql;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.bh"); ++dem;
    values[dem]= -thetaa;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.bh_p"); ++dem;
    values[dem]= deltab*s.ip;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.cb"); ++dem;
    values[dem]= -1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.tau"); ++dem;
    values[dem]= -1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.n_p"); ++dem;
    values[dem]= (Rn-1)*s.ip;
    ++dem;
    
    //g7
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.n"); ++dem;
    values[dem]= 1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.m"); ++dem;
    values[dem]= -varphi;
    ++dem;

    //g8
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.n"); ++dem;
    values[dem]= 1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.bh"); ++dem;
    values[dem]= (1-vec_kappa[t]);
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.m"); ++dem;
    values[dem]= -1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.muc"); ++dem;
    values[dem]= 1;
    ++dem;

    
    //g9
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.bh"); ++dem;
    values[dem]= 1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.bh_p"); ++dem;
    values[dem]= -deltab*s.ip;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.ip"); ++dem;
    values[dem]= -deltab*s.bh_p;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.s"); ++dem;
    values[dem]= -1;
    ++dem;
    
    
    //Households
    
    //g10
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.lama"); ++dem;
    values[dem]= -1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.etaz"); ++dem;
    values[dem]= -2*pow(fmax(s.etaz,0),1);
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.ch"); ++dem;
    values[dem]= -1/s.ch/s.ch;
    ++dem;
    
    //g11
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.lamb"); ++dem;
    values[dem]= s.ch;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.ch"); ++dem;
    values[dem]= s.lamb;
    ++dem;
    
    //g12
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.lama"); ++dem;
    values[dem]= 1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.Rm"); ++dem;
    values[dem]= -beta_h*s.lamb_f*s.ip_f;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.lamb_f"); ++dem;
    values[dem]= -beta_h*s.Rm*s.ip_f;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.ip_f"); ++dem;
    values[dem]= -beta_h*s.Rm*s.lamb_f;
    ++dem;
    
    //g13
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.ql"); ++dem;
    values[dem]= s.lamb;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.lamb"); ++dem;
    values[dem]= s.ql;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.ql_f"); ++dem;
    values[dem]= -beta_h*deltab*s.lamb_f*s.ip_f;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.lamb_f"); ++dem;
    values[dem]= -beta_h*(deltab+ deltab*s.ql_f)*s.ip_f;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.ip_f"); ++dem;
    values[dem]= -beta_h*(deltab+ deltab*s.ql_f)*s.lamb_f;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.etab"); ++dem;
    values[dem]= -2*pe*pow(fmax(s.etab,0),1);
    ++dem;
    
    //g14
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.lamb"); ++dem;
    values[dem]= 1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.lamb_f"); ++dem;
    values[dem]= -beta_h*(1-deltak);
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.pm_f"); ++dem;
    values[dem]= -beta_h*alphaa*s.lama_f*s.y_f/s.k;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.lama_f"); ++dem;
    values[dem]= -beta_h*alphaa*s.pm_f*s.y_f/s.k;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.k"); ++dem;
    values[dem]= -beta_h*alphaa*s.pm_f*s.lama_f*s.y_f*(-1)/s.k/s.k;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.y_f"); ++dem;
    values[dem]= -beta_h*alphaa*s.pm_f*s.lama_f/s.k;
    ++dem;

    
    //g15
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.Rm_p"); ++dem;
    values[dem]= s.m_p*s.ip;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.m_p"); ++dem;
    values[dem]= s.Rm_p*s.ip;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.ip"); ++dem;
    values[dem]= s.Rm_p*s.m_p- deltab*s.bh_p;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.ql"); ++dem;
    values[dem]= s.s;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.s"); ++dem;
    values[dem]= s.ql;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.tau"); ++dem;
    values[dem]= 0;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.ch"); ++dem;
    values[dem]= -1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.i"); ++dem;
    values[dem]= -1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.bh_p"); ++dem;
    values[dem]= -deltab*s.ip;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.etaz"); ++dem;
    values[dem]= 0;
    if (s.etaz < 0){
        values[dem]= -2*s.etaz;
    }
    ++dem;
    
    
    //g16
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.bh"); ++dem;
    values[dem]= -1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.etab"); ++dem;
    values[dem]= 1;
    ++dem;

    
    //g17
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.pie"); ++dem;
    values[dem]= -iota*2*s.pie+ iota;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.lama_f"); ++dem;
    values[dem]= iota*beta_h/s.lama*(s.pie_f-1)*s.pie_f*s.y_f/s.y;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.lama"); ++dem;
    values[dem]= iota*beta_h*s.lama_f*(-1)/s.lama/s.lama*(s.pie_f-1)*s.pie_f*s.y_f/s.y;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.pie_f"); ++dem;
    values[dem]= iota*beta_h*s.lama_f/s.lama*(2*s.pie_f-1)*s.y_f/s.y;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.y_f"); ++dem;
    values[dem]= iota*beta_h*s.lama_f/s.lama*(s.pie_f-1)*s.pie_f/s.y;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.y"); ++dem;
    values[dem]= iota*beta_h*s.lama_f/s.lama*(s.pie_f-1)*s.pie_f*s.y_f*(-1)/s.y/s.y;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.pm"); ++dem;
    values[dem]= epsilon;
    ++dem;
    
    //g18
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.y"); ++dem;
    values[dem]= 1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.k_p"); ++dem;
    values[dem]= -alphaa*pow(s.k_p, alphaa-1)*pow(s.l,1-alphaa);
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.l"); ++dem;
    values[dem]= -(1-alphaa)*pow(s.l,-alphaa)*pow(s.k_p,alphaa);
    ++dem;
    
    
    //g19
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.u"); ++dem;
    values[dem]= 1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.pie_f"); ++dem;
    values[dem]= -1/beta_b*phipie*pow(s.pie_f, phipie-1)*pow(s.y/ss[original_map.at("s.y")], phiy);
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.y"); ++dem;
    values[dem]= -1/beta_b*pow(s.pie_f,phipie)*pow(ss[original_map.at("s.y")],-phiy)*phiy*pow(s.y, phiy-1);
    ++dem;
    
    
    //g20
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.y"); ++dem;
    values[dem]= 1- iota/2*(s.pie-1)*(s.pie-1);
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.cb"); ++dem;
    values[dem]= -1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.ch"); ++dem;
    values[dem]= -1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.i"); ++dem;
    values[dem]= -1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.bh"); ++dem;
    values[dem]= -thetaa;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.pie"); ++dem;
    values[dem]= -iota*s.y*(s.pie-1);
    ++dem;
    
    //g21
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.k"); ++dem;
    values[dem]= 1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.k_p"); ++dem;
    values[dem]= -(1-deltak);
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.i"); ++dem;
    values[dem]= -1;
    ++dem;
    
    
    //g22
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.ip"); ++dem;
    values[dem]= s.pie;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.pie"); ++dem;
    values[dem]= s.ip;
    ++dem;
    
    //g23
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.Rf"); ++dem;
    values[dem]= 1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.u"); ++dem;
    if (s.u > Rf_min){
        double texp = exp(s_max*(Rf_min- s.u));
        values[dem]= - ( 1- texp/(1+texp));
        ++dem;
    }
    else {
        double texp = exp(s_max*(s.u- Rf_min));
        values[dem]= - texp/(1+texp);
        ++dem;
    }
    
    //g24
//    iRow[dem]=pg; jCol[dem]=f.at("s.l"); ++dem;
    values[dem]= chi*(nu+1)*pow(s.l, nu);
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.pm"); ++dem;
    values[dem]= -(1-alphaa)*s.y*s.lama;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.y"); ++dem;
    values[dem]= -(1-alphaa)*s.pm*s.lama;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.lama"); ++dem;
    values[dem]= -(1-alphaa)*s.pm*s.y;
    ++dem;
    
    
}

//***********************************************************************************
//***********************************************************************************


void find_jac_struct_for_initial_condition (Index* iRow, Index* jCol, int& dem)
{
    // Starting position
    int pg = (T-1)*var_foc; //Staring of equation
    int px = 0; // Staring variable;
    
    for (int i=0; i<var_foc; ++i){
        iRow[dem]=pg+i; jCol[dem]=px+i; ++dem;
    }
}

void find_jac_values_for_initial_condition (const Number* x, Number* values, int& dem)
{
    for (int i=0; i<var_foc; ++i){
//        iRow[dem]=pg+i; jCol[dem]=px+i; ++dem;
        values[dem]= 1;
        ++dem;
    }
}

//***********************************************************************************
//***********************************************************************************

void find_jac_struct_for_terminal_condition (Index* iRow, Index* jCol, int& dem)
{
    // Starting position
    int pg = T*var_foc; //Staring of equation
    int px = T*var_foc; // Staring variable;
    
    for (int i=0; i<var_foc; ++i){
        iRow[dem]=pg+i; jCol[dem]=px+i; ++dem;
    }
}

void find_jac_values_for_terminal_condition (const Number* x, Number* values, int& dem)
{
    for (int i=0; i<var_foc; ++i){
//        iRow[dem]=pg+i; jCol[dem]=px+i; ++dem;
        values[dem] = 1;
        ++dem;
    }

}

//***********************************************************************************
//***********************************************************************************

void summarize_jac_struct (Index* iRow, Index* jCol, int& dem)
{
    for (int t=1; t<T; ++t){
        find_jacobian_structure_for_time_t(t, iRow, jCol, dem);
    }
    find_jac_struct_for_initial_condition(iRow, jCol, dem);
    find_jac_struct_for_terminal_condition(iRow, jCol, dem);
}

void summarize_jac_values (const std::vector<double>& ss, const Number* x, Number* values, int& dem)
{
    for (int t=1; t<T; ++t){
        find_jacobian_values_for_time_t(t, ss, x, values, dem);
    }
    find_jac_values_for_initial_condition(x, values, dem);
    find_jac_values_for_terminal_condition(x, values, dem);
}






